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Abstract 

In  recent  years,  there  has  been  a  revival  of  interest  in  the  Schwarz 
alternating  method.  Other  domain  decomposition  algorithms,  in  par- 
ticular the  so  called  iterative  substructuring  methods,  have  also  been 
developed  to  solve  elliptic  finite  element  problems.  In  this  paper,  we 
present  an  additive  variant  of  the  Schwarz  method  for  elliptic  problems, 
which  shows  great  promise  for  parallel  computers.  By  using  techniques 
previously  developed  for  iterative  substructuring  methods,  we  are  able 
to  show  that  this  method  converges  quite  rapidly  even  when  the  re- 
gion is  divided  into  many  subregions.  We  note  that,  as  is  the  case  with 
other  fast  iterative  methods  for  many  subregions,  a  mechanism  for  the 
global  transportation  of  information  is  necessary  in  order  to  obtain  fast 
convergence. 
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1      Introduction 

In  this  paper,  we  will  consider  linear,  self  adjoint,  elliptic  problems  dis- 
cretized  by  finite  element  methods.  Our  purpose  is  to  develop  and  study 
an  additive  variant  of  the  classical  Schwarz  procedure,  see  [16]  ,  which  con- 
verges quite  rapidly  even  for  a  large  number  of  subregions  and  which  permits 
the  use  of  conjugate  gradient  acceleration.  Adopting  a  term  of  structural 
engineering,  we  will  refer  to  certain  non-overlapping  subregions  as  substruc- 
tures. 

For  parallel  computing  systems  with  many  processors  and  in  the  study 
of  speed-up  of  algorithms  on  such  systems,  the  case  of  many  substructures  is 
clearly  of  particular  interest.  The  subproblems  can  then  be  handled  conve- 
niently by  separate  processors  .  If  the  work  of  the  linear  equation  solver  in 
use  grows  faster  than  linearly  with  the  number  of  degrees  of  freedom,  then 
fast  domain  decomposition  methods  can  often  also  offer  real  benefits,  even 
if  only  one  processor  is  used.  Sometimes  we  can  also  use  subregions  with 
simple  geometries  to  reduce  the  overall  effort  of  solving  the  problem. 

In  our  previous  work  on  iterative  substructuring  methods,  we  have  es- 
tablished that  a  mechanism  for  the  global  transportation  of  information  is 
required  in  order  to  obtain  a  rate  of  convergence  that  varies  only  very  slowly 
with  the  number  of  degrees  of  freedom  and  subregions  ;  see  Widlund  [20], 
[19]  .  In  those  papers,  we  have  shown  that  if  we  only  have  next  neighbors 
communication,  the  minimum  number  of  iterations  required  grows  at  least 
as  fast  as  N^^^.  where  N  is  the  number  of  substructures.  The  same  argument 
applies  almost  verbatim  here  and  it  wiU  therefore  not  be  repeated.  In  our 
new  algorithm  for  problems  in  the  plane,  we  use  one  of  the  devices  which  has 
been  successful  in  the  development  of  iterative  substructuring  methods.  Our 
main  result  is  somewhat  stronger  than  those  discussed  in  [19].  The  number 
of  iterations  required,  to  decrease  the  energy  norm  of  the  error  by  a  fixed 
factor,  is  proportional  to  \/{l  +  \og{H/h)).  Here  H  and  h  are  the  diame- 
ters of  a  typical  substructure  and  an  element,  respectively.  The  hierarchical 
bases  method  due  to  Yserentant  [21]  and  an  iterative  substructuring  method 
due  to  Bramble,  Pasciak  and  Schatz  [4]  are  among  the  methods  considered 
in  Widlund  [19]  .  The  bounds  for  these  are  proportional  to  (1  +  log{H/h)). 
We  note  that,  to  our  knowledge,  a  generalization  of  Yserentant's  algorithm 
to  three  dimensional  problems  with  fast  convergence  has  not  been  found, 
while  efficient  iterative  substructuring  Igorithms  have  been  developed  for 
that  case;  cf.  Bramble.  Pasciak  and  Schatz  [5]  and  Dryja  [7]. 

For  self-adjoint  elliptic  problems,  the  finite  element  solution  is  the  pro- 
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jection  of  the  exact  solution  onto  a  finite-dimensional  space.  Different  vari- 
ants of  Schwarz'  algorithm,  in  particular  the  additive  form  considered  in 
this  paper,  can  also  conveniently  be  described  in  terms  of  projections  onto 
subspaces  related  to  overlapping  subregions  which  cover  the  region.  These 
projections  are  orthogonal  with  respect  to  the  symmetric  bihnear  form  as- 
sociated with  the  elliptic  operator  under  consideration.  In  addition,  we  are 
using  a  special  subspace  which  gives  rise  to  a  coarse  finite  element  model 
and  which  provides  the  global  transfer  of  information  that  is  required. 

Our  point  of  departure  has  been  a  recent  paper  by  P.-L.  Lions  [10]  in 
which  a  variational  frame  work  for  the  classical,  multiplicative  Schwarz' 
method  is  developed  for  continuous  elliptic  problems.  Since  it  is  not  nec- 
essary to  rely  on  a  maximum  principle,  the  Schwarz  method  can  be  used 
for  a  number  of  important  elhptic  problems,  such  as  the  equations  of  linear 
elasticity.  Stokes'  equations  and  the  biharmonic  equation.  The  analysis  of 
finite  element  models  is  similarly  made  much  simpler.  In  his  paper.  Lions 
gives  a  number  of  references  to  older  results.  Thus,  more  than  50  years 
ago,  Sobolev  [17]  introduced  a  variational  frame  work  in  a  study  of  Schwarz' 
method  and  the  equations  of  linear  elasticity.  Early  work  was  also  done  by 
Babuska  [1],  [2]  and  by  Morgenstern  [13]:  see  Lions  [10]. 

The  additive  variant  of  the  algorithm  has  been  discovered  independently 
by  us,  but  its  origin  is  not  clear.  We  are  quite  interested  in  finding  out  about 
any  previous  work  of  this  nature.  We  have  found  an  algorithm  quite  similar 
to  ours,  in  a  paper  by  Nepomnyashchikh  [14]  ,  on  iterative  substructuring 
methods  on  regions  subdivided  into  a  fixed  number  of  subregions.  We  also 
note  that  we  have  discovered  that  the  so  called  FAC  and  AFAC  algorithms, 
which  are  iterative  refinement  algorithms  developed  by  McCormick  and  his 
coworkers,  see  [11]  ,[12],  [9],  can  be  described  in  terms  of  projections  and 
multiplicative  and  additive  algorithms,  respectively.  See  Widlund  [18]  for  a 
further  discussion. 

In  Section  2,  we  introduce  our  algorithm  and  prove  that  it  converges 
rapidly  for  Poisson's  equation  discretized  by  piecewise  linear,  continuous 
finite  elements. 


2      The  Additive  Algorithm  and  an  Analysis  of  a 
Model  Problem 

In  this  paper,  we  consider  linear,  self  adjoint,  elliptic  problems  discretized 
by  finite  element  methods  on  Lipschitz  regions.  In  this  section,  we  assume 
that  fi  is  a  plane  region,  that  the  differential  operator  is  the  Laplacian  and 
that  we  use  continuous,  piecewise  linear  finite  elements.  The  continuous  and 
discrete  problems  are  of  the  form 

a{u,v)  =  /(f),  V  r  e    V  , 

and 

a[nh,Vh]  =  f(vh).^  ^  Vh  6    V\  (1) 

respectively.  The  bilinear  form  is  defined  by 


a{u,v)  —    /    Vu-Vvdx. 
Jq 

This  bihnear  form  defines  a  semi-norm  |u|//i  =  (a(u,u))^/^  in  H^{Q  ).  We 
assume  that  dQ.  the  boundary  of  0  .  is  the  union  of  two  non-overlapping 
sets,  Td  and  Fa',  on  which  zero  Dirichlet  and  arbitrary  Neumann  bound- 
ary conditions  are  given.  We  assume  that  T p  is  non  empty,  to  avoid  an 
unnecessary  complication  with  a  singular  problem.  All  elements  of  V  and 
its  subspace  V^  vanish  on  T d-  ^'e  note  that  we  can  always  reduce  an  inho- 
mogenous  Dirichlet  problem  to  the  case  considered  here  at  the  expense  of 
solving  one  problem  on  each  of  the  subregions  which  are  next  to  To  ■  The 
subregions  are  introduced  below. 

We  introduce  the  triangulation  of  J7  in  the  following  way.  We  first  divide 
the  region  into  non-overlapping  substructures  fi,,  z  =  1,-  ■  ■  ,N  .  To  simplify 
the  description  and  analysis  of  our  algorithm  ,  we  confine  our  study  to  tri- 
angular substructures.  In  such  a  case,  the  original  region  must  of  course  be 
a  polygon  .  We  also  adopt  the  common  assumption  in  finite  element  theory 
that  all  substructures  are  shape  regular  in  the  sense  that  the  diameter  Hi, 
divided  by  the  radius  of  the  largest  inscribed  circle  in  Q,  ,  are  uniformly 
bounded.  All  the  substructures  Q,  are  subdivided  into  elements.  The  ele- 
ments are  shape  regular  in  the  same  sense  as  above.  Since  the  Schwarz-type 
domain  decomposition  algorithms  use  overlapping  subregions,  we  extend 
each  substructure  to  a  larger  region  Q'-  .  We  assume  that  the  distance 
between  the  boundaries  dfl,  and  dQ[  is  bounded  from  below  by  a  fixed 
fraction  of  H,,  and  that  dil[    does  not  cut  through  any  elment.   We  make 


the  same  construction  for  the  substructures  that  meet  the  boundary  except 
that  we  cut  off  the  part  of  Q[  that  is  outside  of  fi  . 

We  remark  that  some  of  the  more  intricate  issues  in  the  analysis  of  the 
Schwarz  methods  arise  when  the  boundaries  of  the  different  subdomains  Q[ 
intersect  at  one  or  several  points;  cf.  the  discussion  in  Lions  [10].  This 
happens  for  example  if  the  region  is  L-shaped  and  is  partitioned  into  two 
overlapping  rectangles.  Tools  powerful  enough  to  treat  problems  of  this  kind 
are  known  to  us  .  Sometimes  the  resulting  bounds  involve  an  extra  factor 
(1  +  \og{H/h))  ;  see  further  Widlund  [18]  .  In  this  paper  ,  we  will  confine 
our  study  to  the  simpler  case  . 

Our  finite  element  space  is  represented  as  the  sum  of  N  +  1  subspaces 

The  first  subspace  \'q  .  which  we  also  call  V^,  is  special  .  It  is  the  space 
of  continuous,  piecewise  linear  functions  on  the  coarse  mesh  defined  by  the 
substructures  0,  .  When  convenient,  we  can  view  any  element  of  V  as  a 
coarse  mesli  interpolant  ///i>.  of  some  element  v^  €  V"  .  The  computation 
of  the  projection  of  an  arbitrary  function  onto  this  subspace  involves  the 
solution  of  a  standard  finite  element  linear  system  of  algebraic  equations 
which  is  on  the  order  of  N.  This  coarse,  global  approximation  of  the  ellip- 
tic equation  can  often  practically  be  handled  by  a  direct  method  such  as 
Gaussian  elimination. 

The  subspaces,  \\^  ,  which  are  related  to  interior  substructures  ,  are 
defined  as  1'''  f]  Hl{n[)  .  The  space  Hq  is,  as  usual,  the  subspace  of  H^ 
functions  with  zero  trace.  It  is  well  known  that  this  space  can  be  extended 
continuously  by  zero.  After  extending  Hq(^'j)  by  zero  ,  we  can  regard 
V-  as  a  subspace  of  \  .  We  construct  subspaces  related  to  the  boundary 
substructures  similarly.  The  boundary  conditions  of  the  original  problem 
are  inherited  at  any  node  of  fi'  which  falls  on  dfl.  A  nodal  value  is  thus 
constrained  to  be  zero  if  the  node  belongs  to  Tp  ,  while  it  is  free  at  any  node 
on  T^'  .  At  all  nodes  of  dQ\  which  are  in  the  interior  of  fl  we  impose  zero 
Dirichlet  conditions. 

The  projection  Pyh  =  P,  is  defined,  for  all  of  V'^,  by  the  unique  element 
of  T^'',  which  satisfies 

aiP,Vh,4>h)  =  aivh,4>h),  '^4>h  e  Vt.  (2) 

Lions  [10]  has  shown  that  the  error  propagation  operator  of  the  standard 
multiplicative  variant  of  the  Schwarz  method  can  be  written  as 

(/-Pl)(/-P2), 


for  the  case  of  two  subregions.  We  can  view  that  algorithm  as  an  iterative 
method  for  solving 

with  an  appropriate  right  hand  side  g^,  ■  We  note  that  the  operator  is  a  poly- 
nomial of  degree  two,  and  thus  is  not  ideal  for  parallel  computing,  since  two 
sequential  steps  are  involved  .  The  additive  form  of  the  Schwarz  algorithm 
can  similarly  be  regarded  as  an  iterative  method  for  solving  the  equation 

Pu^,  =  (Po  +  Pi  +  ■■■■  +  PN)^h  =  g'k  .  (3) 

W'e  must  assure  that  this  equation  has  the  same  solution  as  equation  1  ,  i.e. 
we  must  use  the  correct  right  hand  side.  Since  by  equation  1,  we  have 

a[uh,(ph)  -  f[4)h)  , 

we  can  construct  the  right-hand  side  5^  by  solving  equation  2  for  all  values 
of  ;'  and  adding  the  resulting  vectors  .  We  can  also  apply  the  operator  P 
of  equation  3  to  any  given  element  of  V''  at  the  same  expense  by  applying 
each  projection  P,  once  and  adding  the  vectors  .  We  solve  equation  3  by  the 
standard  conjugate  gradient  method.  Much  of  the  work,  in  particular  that 
involving  the  individual  projections,  can  be  carried  out  in  parallel  . 

It  is  well  known  that  the  number  of  steps  required  to  decrease  an  appro- 
priate norm  of  the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor 
is  proportional  to  ^/K  ,  where  k  is  the  condition  number  of  P  ;  see  .  e.g. 
Golub  and  Van  Loan  [8]  .  We  therefore  need  to  establish  that  the  opera- 
tor P  of  equation  3  is  not  only  invertible  but  that  satisfactory  upper  and 
lower  bounds  on  its  eigenvalues  can  be  obtained.  We  will  prove  the  following 
theorem. 

Theorem  1    n{P)    is  bounded  by  const  (1  -f  log{H/h)). 

Here,  and  elsewhere  in  this  paper,  the  constants  are  independent  of  the 
diameters  of  of  the  substructures  and  the  individual  elements.  They  vary 
with  the  shape  of  the  elements  and  substructures.  Thus  in  two  dimensions, 
they  depend  on  the  minimum  angle.  As  will  become  apparent  in  the  proof, 
the  constant  in  the  theorem  also  grows  if  the  overlap  of  the  subregions 
shrinks. 

The  upper  bound  on  the  spectrum  is  obtained  by  bounding 

aiPvh-,Vh)  =  a{PoVh,Vh)  +  a{PiVh,Vh)  +  ■  ■  ■  +  a{P^'Vh,Vh)  , 


from  above  in  terms  of  a(t'/,,t'/i)  .  We  first  use  Schwarz'  inequality  and 
the  fact  that  Pq  is  a  projection  to  prove  that  the  first  term  is  bounded  by 
aivh,Vk)  .  Similarly  ,  it  is  easy  to  show  that  the  spectrum  of  P  is  bounded 
from  above  by  (A'  +  1)  .  Our  goal,  however,  is  to  establish  a  bound  on  the 
condition  number  that  grows  quite  slowly,  and  we  therefore  need  a  better 
upper  bound  . 

The  sum  of  the  other  terms,  which  are  of  the  form 

a{P,Vh,Vh)  =  a(/',r/,,P,t;/i)  , 

can  also  be  estimated  by  const  a{vh.Vh)  ■  To  see  this  .  we  will  use  the  fact 
that  PiVh  vanishes  outside  n[  and  that  there  is  a  uniform  bound  on  the 
number  of  terms  that  are  different  from  zero  at  any  given  point  .  It  is  also 
easy  to  see,  from  the  definition  of  the  \']^  ,  that  P,Vk  depends  only  on  the 
values  of  r/i    in  Q[  .  Therefore 

P,Vh  -  P,ii',Vh)  , 

where  V'l  is  a  smooth  cut-off  function.  It  has  values  between  zero  and  one,  is 
equal  to  one  on  fi'  and  vanishes  outside  a  neighborhood  Q"  of  fij  .  We  can 
choose  this  new  set  so  that  the  distance  between  its  boundary  and  that  of 
(7'  is  on  the  order  of  H,  ,  the  diameter  of  fi,  .  We  modify  this  construction 
for  boundary  substructures  by  cutting  off  the  part  that  is  outside  fi  .  It  is 
then  easy  to  construct  V'i  so  that  its  gradient  is  bounded  by  const  j Hi  . 
By  using  the  formula  for  the  derivative  of  a  product,  we  then  find  that  the 
solution  of  equation  2  satisfies 

(P.t'A.P.n)Hi(fi)  =  const  {\vk\''^^^^,,^  +  H-^  IK'/vlli,(n;'))  •  (4) 

We  wish  to  remove  the  L2  -  norm  term  from  the  estimate  4  and  we 
have  to  consider  two  cases.  If  5fi"  does  not  intersect  T/j,  we  first  modify  u/, 
locally  by  adding  an  arbitrary  constant  .  This  will  leave  the  right  hand  side 
unchanged,  since  as  is  easily  shown  by  using  the  definition  of  the  projection 
that  Piconst  =  0.  The  semi-norm  in  the  right  hand  side  of  4  also  remains 
unchanged.  Poincare's  inequahty,  see  e.g.  Necas  [15]  ,  can  now  be  used 
to  estimate  the  second  term  of  the  right  hand  side  by  the  first.  If,  on  the 
other  hand,  dVl'l  intersects  Yd  ,  we  cannot  shift  by  a  constant  and  remain 
in  the  space  V  .  Instead,  we  use  an  inequality  due  to  Friedrichs,  given  as 
Theorem  1.9  in  Necas  [15]  ,  to  estimate  the  L2{Cl")  -  term  .  This  inequality 
essentially  states  that,  for  any  region  (l  of  diameter  one, 

\Ml^^n)<  const  {\u\]j,^^^  +  \\u\\l^^^)). 


Here  7  is  a  subset  of  positive  measure  of  dQ  .  If  we  scale  the  variables,  we 
find  that  the  constant  in  the  estimate  decreases  in  proportion  to  Hf  ,  if  the 
measure  of  dCl[  f]  Td  remains  bounded  from  below  by  that  of  dfl[  .  This 
latter  condition  causes  no  difficulty,  since,  if  necessary,  we  can  always  make 
fi"  somewhat  larger.  The  factor  Hf  is  exactly  what  is  required  to  control 
the  X2(n")  -  term  .   Our  proof  of  the  upper  bound  for  P  is  completed  by 

noting  that 

N 

since  the  fi"  provides  a  finite  cover  of  fi  . 

We  note  that  we  have  essentially  used  a  quotient  space  argument,  which 
also  plays  an  important  role  in  the  error  analysis  for  finite  element  methods; 
cf.  Ciarlet  [6]  .  Similar  arguments  have  been  used  extensively  in  the  theory 
of  iterative  substructuring  methods;  see  e.g.  Widlund  [20],  [19]  . 

We  use  a  lemma  in  Lions  [10]  to  obtain  a  lower  bound.  Since  the  proof 
of  his  result  is  quite  short,  we  include  it  in  this  paper. 

Leinma  1  (Lions)  Lei  u^,  =  J2i=o'^h.i  ■,  U'here  u/,,,  6  V',''  ,  be  a  partition  of 
an  element  of\''^'   and  assume  further  that  XIi=o  l"/i,i|j/i  ^  Colu/ij^i  ,  Vu/i  G 

-0 


V^  .    Then  Xm:n{P)  >  Co"' 


Proof:   By  elementary  properties  of  symmetric  projections  and  the  rep- 
resentation of  vt,   as  a  sum,  we  find  that 

A'  A'  A' 

t=0  1=0  t=0 

Therefore, 

t=0  1=0 

By  the  assumption  of  the  lemma 

A'  N 

Whln^  <  ClJ2\P'''h\H^  =  C2^(P,UA,ti/.)Hi  =  C^(Puh,Uh)„.  , 

t=0  t=0 

and  the  lemma  is  established. 

What  remains  is  to  find  a  suitable  decomposition  and  to  establish  the 
bound  required  by  Lions'  lemma.  We  first  write 

Uh  =  lu^h  +  {uh  -  Ihuh)  =  Ihuh  +  Wh  , 


where  ///  is  the  interpolation  operator  related  to  the  space  V''-'^  =  Vq  ,  which 
was  introduced  above.  We  use  the  following  lemma,  which  plays  a  similar 
role  in  the  theory  for  iterative  substructuring  algorithms.  This  result  ,  which 
dates  back  to  at  least  19G6  ,  is  proven  in  a  number  of  papers;  see  e.g.  Bramble 
[3]    ,  Bramble,  Pasciak  and  Schatz  [4]  or  Yserentant  [21]    . 

Lemma  2 


It  is  easy  to  see  that  |///i//^,p  j  q,    can  be  estimated  by  a  sum  of  differ- 

ences  of  functions  values  of  Uh  at  the  vertices  of  Cl[  .  As  shown  in  Widlund 
[20],  [19],  we  can  use  the  lemma  together  with  the  same  quotient  space 
argument  .  as  used  to  derive  our  upper  bound  ,  to  establish  that 

|///"/.r;^,,0;)  ^  <^onst{l  +  logiH/h))\uh\]^^^^,^^  .  (5) 

The  other  terms  in  the  decomposition  are  defined  by  u/j,,  =  //i(^,i(.'/i),7  = 
0,  •  ■  • ,  jV  .  Here  the  6,  define  a  partion  of  unity  and  belong  to  Cq^{Q[)  .  We 
can  arrange  so  that  V^,  is  bounded  by  const  /H,  .  By  using  the  linearity  of 
Ifi  ,  we  can  easily  show  that  we  obtain  a  correct  partitioning  of  u/;  .  In  order 
to  estimate  the  semi-norm  of  Uh,i  ,  we  partition  fl[  into  elements.  For  each 
element  K.  we  obtain 

Here  6,  is  the  average  value  of  ^,  over  K.  It  is  easy  to  see,  by  using  an  inverse 

inequality,  that 

\hiie.-r,)wh)\H^iiq  <  const  h-'\\h{{e,  -F,)wh)\\l,(K)  ■ 

We  can  now  use  the  fact  that  on  K,  6,  differs  from  its  average  by  at  most 
const  h/H,  ,  and,  after  summing  over  all  elements  of  Cl'-  ,  we  arrive  at  an 
inequality  similar  to  4: 

At  the  expense  of  introducing  a  factor  (1  +  log{H/h))  ,  we  can  estimate  the 
right  hand  side  of  6  by  the  same  expression  with  Uf,  +  const  instead  of  w^  . 
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Here  the  constant  term  is  arbitrary.  We  first  note  that  w^  =  [I  —  lH)uh  does 
not  change  if  we  add  a  constant  to  i//,  .  By  the  inequality  5,  we  can  estimate 
the  first  term  of  the  right  hand  side  of  6  as  required.  The  second  term  can 
be  estimated  directly  by  ||u-;,||'^^  q/    ,  since  the  area  of  fij  is  on  the  order 

of  H}  .  This  quantity  can  easily  be  estimated  by  \\uh  +  const  W'l^.n,/.  ■  We 

now  use  Lemma  2  ,  Poincare's  inequality  and  a  finite  covering  argument  to 
complete  the  proof  of  our  Theorem. 
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